This script is to compare model fits between 1) the old model (unpooled K), 2) the pooled-K model we developed, with pooled and logged K, and 3) a non-centered, pooled-K model.
Summary: the new pooled-logK model is still really struggling for two sites: BIGC and CARI. ChatGPT and other resources suggested a non-centered model structure to avoid ‘funnelling’ for sigma and mu: instead of drawing each logU or logK directly from a normal distribution with unknown σ, the non-centered version draws a standard-normal z first, then scales and shifts it by μ and σ.
For BIGC, the new data is part of the problem: the earlier, unpooled-K model is also struggling. The non-centered model works well. Solution: find days that might be causing the problems, remove from the dataset, try all models again.
For CARI, the original unpooled-K model runs well, but not the initial pooled-K model. The uncentered pooled-K model works well.
Parameter values are similar across all three models, which is a good sign.
For the other four sites, CUPE, PRIN, SYCA, WLOU:
The original pooled-K model does very well for CUPE and SYCA, and decently well for PRIN and WLOU. More iterations will probably bring all r-hats to 1.
The non-centered pooled-K model is supposed to improve model performance across the board, but it does not. As noted, it worked really well for BIGC and CARI, decently well but more slowly (vs the other pooled-K model) for most of the others (CUPE, PRIN, SYCA), and poorly for WLOU. Again, parameter estimates are similar across all three models (unpooled, pooled logK, uncentered pooled logK), so not much to choose from there.
BIGC probably has problem days, removing them will hopefully improve the modeling. CARI does well w. the non-centered model and struggles with our initial pooled-K (‘centered’) model. WLOU does well w. the ‘centered’ model and struggles with the non-centered model. Thoughts on how to proceed?
Fits, chains follow
# Load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstanarm)
library(here)
library(plotly)
library(sf)
library(MetBrewer)
library(PNWColors)
library(ggpubr)
library(cowplot)
library(ragg)
library(latex2exp)########## Load model data _____________________________________________________
# The same model data should work for both pooled and unpooled models
bigc.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/bigc.data.rds"))
cari.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cari.data.rds"))
cupe.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cupe.data.rds"))
prin.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/prin.data.rds"))
syca.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/syca.data.rds"))
wlou.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/wlou.data.rds"))########## Load unpooled model fit ______________________________________________________
bigc.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/bigc.fit.rds"))Print unpooled parameter fit
print(bigc.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 2049 1.00
sigma_U 0.07 0.01 0.03 0.02 0.05 0.07 0.09 0.15 15 1.28
b0 -7.30 0.13 0.94 -9.29 -7.92 -7.25 -6.65 -5.55 53 1.07
b1 0.41 0.01 0.08 0.26 0.36 0.41 0.47 0.59 53 1.07
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:23:46 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(bigc.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)
########## Load pooled logK model fit ______________________________________________________
bigc.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/bigc.fit.rds"))Print pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(bigc.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 746 1.00
sigma_U 0.08 0.01 0.03 0.02 0.05 0.07 0.10 0.14 15 1.30
sigma_logK 0.44 0.01 0.07 0.30 0.39 0.44 0.49 0.59 74 1.04
logK_mean 0.83 0.01 0.09 0.65 0.78 0.83 0.89 1.00 75 1.04
b0 -7.39 0.36 1.29 -10.45 -8.06 -7.20 -6.45 -5.33 13 1.32
b1 0.42 0.03 0.11 0.23 0.34 0.40 0.48 0.69 13 1.32
Samples were drawn using NUTS(diag_e) at Tue Oct 28 09:53:24 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(bigc.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)
##### Examine distributions of struggling params +
# # gets fit into a dataframe - easier for the bayesplotting
# posteriorMA <- as_draws_matrix(posterior)
# sigma <- posteriorMA[, "sigma"]
#
# # Find all conc_tilde columns
# conc_cols <- grep("^conc_tilde\\[", colnames(posteriorMA), value = TRUE)
#
# cors <- sapply(conc_cols, function(nm) {
# x <- posteriorMA[, nm]
# if (is.numeric(x) && all(!is.na(x))) {
# cor(x, sigma)
# } else {
# NA
# }
# })
#
# cors <- cors[!is.na(cors)]
# top <- names(sort(abs(cors), decreasing = TRUE))[1]
#
# cat("Most correlated:", top, " (r =", round(cors[top], 3), ")\n")
#
#
# bayesplot::mcmc_scatter(as.array(bigc.fit_pooledK.log),
# pars = c(top, "sigma"),
# transform = list(sigma = "log"),
# alpha = 0.3)
#
#
#
# scatter_sigma_pool <- mcmc_scatter(
# posterior,
# pars = c("sigma", "conc_tilde[10,20]"),
# transform = list(sigma = "log"),
# size = 1
# )
#
# scatter_sigma_pool
# A second, non-centered pooled model worked better: - Model finished in 7 min 45 sec
########## Load pooled logK model fit ______________________________________________________
bigc.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/bigc.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(bigc.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.06 0.00 0.00 0.06 0.06 0.06 0.06 0.06 2320 1.00
sigma_U 0.04 0.00 0.03 0.00 0.02 0.04 0.06 0.12 765 1.00
sigma_logK 0.44 0.00 0.08 0.30 0.38 0.43 0.49 0.59 266 1.02
logK_mean 0.82 0.00 0.09 0.64 0.76 0.82 0.88 0.98 399 1.01
b0 -6.99 0.02 0.97 -8.85 -7.61 -6.99 -6.32 -5.16 1682 1.00
b1 0.38 0.00 0.09 0.22 0.32 0.38 0.44 0.55 1674 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 23:35:13 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(bigc.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - all chains converged well
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 4, size = 2)########## Load unpooled model fit ______________________________________________________
cari.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/cari.fit.rds"))Print unpooled parameter fit
print(cari.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.13 0.00 0.00 0.12 0.13 0.13 0.13 0.13 4363 1
sigma_U 0.40 0.00 0.03 0.34 0.37 0.40 0.42 0.47 579 1
b0 -1.95 0.02 0.53 -3.00 -2.27 -1.94 -1.59 -0.94 1122 1
b1 0.12 0.00 0.05 0.03 0.09 0.12 0.15 0.21 1111 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:37:16 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cari.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines) - not needed when chains look good
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)########## Load pooled logK model fit ______________________________________________________
cari.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cari.fit.rds"))Print pooled logK parameter fit
print(cari.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.12 0.00 0.00 0.12 0.12 0.12 0.12 0.13 1835 1.00
sigma_U 0.38 0.00 0.03 0.32 0.36 0.38 0.40 0.44 403 1.02
sigma_logK 0.07 0.02 0.05 0.01 0.02 0.06 0.10 0.17 4 2.01
logK_mean -1.00 0.01 0.05 -1.10 -1.04 -1.00 -0.97 -0.90 50 1.08
b0 -2.33 0.02 0.50 -3.31 -2.67 -2.31 -1.98 -1.40 927 1.00
b1 0.15 0.00 0.04 0.07 0.12 0.15 0.18 0.24 856 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 10:10:18 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cari.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)
A second, non-centered pooled model composition worked better:
########## Load pooled logK model fit ______________________________________________________
cari.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cari.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(cari.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.12 0.00 0.00 0.12 0.12 0.12 0.12 0.13 2890 1
sigma_U 0.38 0.00 0.03 0.32 0.36 0.38 0.40 0.45 926 1
sigma_logK 0.06 0.00 0.04 0.00 0.02 0.05 0.08 0.17 780 1
logK_mean -1.00 0.00 0.05 -1.10 -1.04 -1.00 -0.97 -0.91 1202 1
b0 -2.32 0.01 0.50 -3.32 -2.65 -2.32 -1.97 -1.36 1201 1
b1 0.15 0.00 0.04 0.07 0.12 0.15 0.18 0.24 1184 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 23:52:03 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cari.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - all chains converged well
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)########## Load unpooled model fit ______________________________________________________
cupe.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/cupe.fit.rds"))Print unpooled parameter fit
print(cupe.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.08 0.00 0.00 0.08 0.08 0.08 0.08 0.08 2730 1
sigma_U 0.30 0.00 0.01 0.27 0.29 0.30 0.30 0.32 2573 1
b0 -3.08 0.01 0.62 -4.29 -3.50 -3.09 -2.67 -1.83 1816 1
b1 0.13 0.00 0.06 0.02 0.09 0.13 0.17 0.24 1821 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 20:51:59 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cupe.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines) - good convergence, did not need
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)########## Load pooled logK model fit ______________________________________________________
cupe.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cupe.fit.rds"))Print pooled logK parameter fit
print(cupe.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.08 0.00 0.00 0.08 0.08 0.08 0.08 0.08 3961 1
sigma_U 0.29 0.00 0.01 0.26 0.28 0.29 0.29 0.31 4110 1
sigma_logK 0.21 0.00 0.01 0.19 0.20 0.21 0.22 0.23 1345 1
logK_mean 1.40 0.00 0.01 1.37 1.39 1.40 1.40 1.42 3143 1
b0 -3.13 0.01 0.61 -4.31 -3.55 -3.12 -2.71 -1.93 3134 1
b1 0.13 0.00 0.06 0.02 0.10 0.13 0.17 0.24 3148 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 10:49:45 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cupe.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - converged fairly well, getting a better look at 1 & 2
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)
A second, non-centered pooled model composition also worked for CUPE, but not quite as well as the original pooled-K - and took a VERY long time, maybe partly due to other draws on computation power/ bandwidth
########## Load pooled logK model fit ______________________________________________________
cupe.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/cupe.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(cupe.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.08 0.00 0.00 0.08 0.08 0.08 0.08 0.08 5206 1.00
sigma_U 0.29 0.00 0.01 0.27 0.28 0.29 0.29 0.31 1015 1.00
sigma_logK 0.21 0.00 0.01 0.18 0.20 0.21 0.22 0.23 1782 1.00
logK_mean 1.40 0.00 0.01 1.37 1.39 1.40 1.40 1.42 2808 1.00
b0 -3.10 0.03 0.61 -4.28 -3.50 -3.10 -2.69 -1.94 466 1.02
b1 0.13 0.00 0.06 0.03 0.09 0.13 0.17 0.24 458 1.02
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:16:26 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(cupe.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - all chains converged well
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 4, size = 2)########## Load unpooled model fit ______________________________________________________
prin.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/prin.fit.rds"))Print unpooled parameter fit
print(prin.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.29 0.00 0.00 0.28 0.29 0.29 0.29 0.30 2841 1
sigma_U 0.49 0.00 0.04 0.42 0.46 0.49 0.51 0.57 976 1
b0 -11.87 0.05 1.69 -15.17 -13.01 -11.83 -10.69 -8.58 980 1
b1 0.98 0.00 0.15 0.68 0.87 0.97 1.08 1.27 988 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:02:43 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(prin.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines) - chains converged, did not need
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)########## Load pooled logK model fit ______________________________________________________
prin.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/prin.fit.rds"))Print pooled logK parameter fit
print(prin.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.28 0.00 0.00 0.28 0.28 0.28 0.29 0.29 3423 1.00
sigma_U 0.48 0.00 0.04 0.41 0.45 0.48 0.50 0.56 721 1.00
sigma_logK 0.44 0.01 0.12 0.24 0.36 0.43 0.51 0.70 80 1.03
logK_mean -0.77 0.01 0.11 -1.02 -0.84 -0.76 -0.69 -0.59 100 1.03
b0 -11.37 0.05 1.62 -14.67 -12.44 -11.35 -10.26 -8.13 1070 1.00
b1 0.93 0.00 0.15 0.64 0.83 0.93 1.03 1.23 1076 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 11:07:32 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(prin.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)
A second, non-centered pooled model composition also worked, but took longer to converge:
########## Load pooled logK model fit ______________________________________________________
prin.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/prin.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(prin.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.28 0.00 0.00 0.28 0.28 0.28 0.29 0.29 2722 1.00
sigma_U 0.47 0.00 0.04 0.40 0.44 0.47 0.50 0.55 609 1.00
sigma_logK 0.39 0.01 0.13 0.12 0.31 0.39 0.48 0.66 152 1.02
logK_mean -0.74 0.01 0.11 -0.98 -0.81 -0.73 -0.66 -0.56 182 1.01
b0 -11.32 0.08 1.66 -14.54 -12.45 -11.33 -10.17 -8.09 434 1.00
b1 0.93 0.01 0.15 0.64 0.82 0.93 1.03 1.21 425 1.00
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:30:56 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(prin.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - all chains converged fairly well with a couple of wobbles in the logK hyperparameters
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)
########## Load unpooled model fit ______________________________________________________
syca.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/syca.fit.rds"))Print unpooled parameter fit
print(syca.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.43 0.00 0.01 0.42 0.43 0.43 0.44 0.45 5439 1
sigma_U 0.75 0.00 0.07 0.63 0.70 0.74 0.79 0.89 2799 1
b0 32.54 0.04 2.91 26.78 30.60 32.60 34.50 38.27 4423 1
b1 -3.77 0.01 0.33 -4.43 -4.00 -3.78 -3.55 -3.12 4358 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:20:27 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(syca.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines) - chains converged, did not need
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "b0", "b1"), highlight = 2, size = 2)########## Load pooled logK model fit ______________________________________________________
syca.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/syca.fit.rds"))Print pooled logK parameter fit
print(syca.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.42 0.00 0.01 0.41 0.42 0.42 0.43 0.43 3119 1.00
sigma_U 0.62 0.00 0.05 0.53 0.58 0.62 0.66 0.73 2148 1.00
sigma_logK 0.18 0.00 0.03 0.12 0.16 0.18 0.20 0.25 465 1.01
logK_mean 2.08 0.00 0.03 2.02 2.06 2.08 2.11 2.15 949 1.00
b0 26.15 0.05 2.44 21.30 24.51 26.17 27.75 30.90 2471 1.00
b1 -3.00 0.01 0.28 -3.54 -3.18 -3.01 -2.81 -2.45 2440 1.00
Samples were drawn using NUTS(diag_e) at Tue Oct 28 11:23:26 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(syca.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - chains converged, did not need
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 4, size = 2)A second, non-centered pooled model composition also worked, but took longer to converge:
########## Load pooled logK model fit ______________________________________________________
syca.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/syca.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(syca.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.42 0.00 0.01 0.41 0.42 0.42 0.43 0.43 4186 1.00
sigma_U 0.62 0.00 0.05 0.53 0.59 0.62 0.66 0.73 571 1.01
sigma_logK 0.18 0.00 0.03 0.12 0.16 0.18 0.20 0.25 884 1.00
logK_mean 2.08 0.00 0.03 2.02 2.06 2.08 2.11 2.15 942 1.00
b0 25.89 0.14 2.43 21.16 24.22 25.95 27.45 30.55 298 1.00
b1 -2.97 0.02 0.28 -3.51 -3.15 -2.98 -2.78 -2.44 311 1.00
Samples were drawn using NUTS(diag_e) at Wed Oct 29 01:57:05 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(syca.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines) - mb a little wobbly in the betas(?)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)
########## Load unpooled model fit ______________________________________________________
wlou.fit_unpooledK <- readRDS(here("N_uptake_NEON/data/model_fits/unpooledK/wlou.fit.rds"))Print unpooled-K parameter fit
print(wlou.fit_unpooledK, pars=c("sigma", "sigma_U", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 9714 1
sigma_U 0.64 0.00 0.04 0.57 0.62 0.64 0.67 0.72 4044 1
b0 -6.74 0.01 0.89 -8.48 -7.34 -6.74 -6.14 -4.98 4338 1
b1 0.43 0.00 0.09 0.25 0.37 0.43 0.49 0.61 4341 1
Samples were drawn using NUTS(diag_e) at Tue Oct 28 21:41:19 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(wlou.fit_unpooledK) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "b0", "b1"))
# highlights one chain (uses points not lines) - chains converged well, did not need
# mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2)########## Load pooled logK model fit ______________________________________________________
wlou.fit_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/wlou.fit.rds"))Print pooled logK parameter fit
print(wlou.fit_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 1582 1.00
sigma_U 0.64 0.00 0.04 0.57 0.61 0.63 0.66 0.71 1269 1.00
sigma_logK 0.42 0.00 0.04 0.33 0.39 0.41 0.45 0.50 237 1.01
logK_mean 1.20 0.00 0.05 1.09 1.16 1.20 1.24 1.30 337 1.00
b0 -6.87 0.03 0.92 -8.71 -7.47 -6.87 -6.24 -5.09 1224 1.00
b1 0.45 0.00 0.10 0.26 0.38 0.45 0.51 0.64 1209 1.00
Samples were drawn using NUTS(diag_e) at Mon Oct 27 17:38:58 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(wlou.fit_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)
A second, non-centered pooled model composition did not work well for WLOU:
########## Load pooled logK model fit ______________________________________________________
wlou.fit2_pooledK.log <- readRDS(here("N_uptake_NEON/data/model_fits/pooledK_logK/wlou.fit2.rds"))Print non-centered pooled logK parameter fit
########## Load pooled logK model fit ______________________________________________________
print(wlou.fit2_pooledK.log, pars=c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))Inference for Stan model: anon_model.
4 chains, each with iter=2500; warmup=1250; thin=1;
post-warmup draws per chain=1250, total post-warmup draws=5000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 0.05 0.00 0.00 0.05 0.05 0.05 0.05 0.05 3066 1.00
sigma_U 0.64 0.00 0.04 0.57 0.61 0.64 0.66 0.71 1549 1.00
sigma_logK 0.48 0.00 0.07 0.35 0.43 0.48 0.53 0.65 330 1.01
logK_mean 1.14 0.00 0.08 0.95 1.09 1.15 1.19 1.27 395 1.01
b0 -6.84 0.03 0.89 -8.63 -7.43 -6.84 -6.23 -5.16 964 1.00
b1 0.44 0.00 0.09 0.27 0.38 0.44 0.51 0.63 969 1.00
Samples were drawn using NUTS(diag_e) at Mon Nov 24 13:50:48 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Examine MCMC chains
posterior <- as.array(wlou.fit2_pooledK.log) # from rstan package
color_scheme_set("viridis")
mcmc_trace(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"))
# highlights one chain (uses points not lines)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 1, size = 2, alpha=0.5)
mcmc_trace_highlight(posterior, pars = c("sigma", "sigma_U", "sigma_logK", "logK_mean", "b0", "b1"), highlight = 2, size = 2, alpha=0.5)